Antiferromagnetism in the Hubbard Model on the Bernal-Stacked Honeycomb Bilayer 
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Using a combination of quantum Monte Carlo simulations, functional renormalization group cal- 
culations and mean-field theory, we study the Hubbard model on the Bernal-stacked honeycomb 
bilayer at half-filling as a model system for bilayer graphene. The free bands consisting of two Fermi 
points with quadratic dispersions lead to a finite density of states at the Fermi level, which triggers 
an antiferromagnetic instability that spontaneously breaks sublattice and spin rotational symmetry 
once local Coulomb repulsions are introduced. Our results reveal an inhomogeneous participation 
of the spin moments in the ordered ground state, with enhanced moments at the threefold coordi- 
nated sites. Furthermore, we find the antiferromagnetic ground state to be robust with respect to 
enhanced interlayer couplings and extended Coulomb interactions. 
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PACS numbers: 71.27. +a,71.10.Fd,71.30.+h,73.21.Ac,75.70.Cn 



There is currently significant interest in understanding 
the electronic properties of bilayer graphene (BLG), in 
particular the ground state at the charge neutrality point. 
Several experimental studies [3-S| hint to the formation 
of a symmetry broken state in BLG, but its actual nature 
remains ambiguous and is at the moment a highly de- 
bated topic. Symmetry breaking in BLG can arise due to 
thermal annealing-induced strain on suspended samples 
as well as external electric fields applied perpendicular to 
the BLG sheets, fn the absence of such external pertur- 
bations, due to the finite density of states at the Fermi 
level in the free band limit, the electronic Coulomb inter- 
action is expected to trigger a genuine electronic instabil- 
ity and drive BLG into a correlated ground state 
sible candidate states that have been suggested 
include an (layered) antiferromagnetic (AF) state, sev- 
eral topological states such as quantum anomalous Hall, 
quantum spin Hall (QSH) or quantum valley Hall states, 
all of which exhibit a finite bulk gap, as well as a gaplcss 
nematic state. While most recent experiments identified 
a finite excitation gap of a few meV emerging in BLG 
at low temperatures 0-11] , the transport data in Ref. Q 
have been interpreted towards the formation of a gapless, 
possibly nematic state. Within the currently inconclusive 
experimental situation, an AF state is considered a prob- 
able ground state [22 - 24 among the (gapfull) candidates 
and thus worth a more detailed examination. Further- 
more, the validity of approximative approaches need to 
be tested against unbiased and numerically exact results. 

Here, we explore the nature of this possible ground 
state by taking screened Coulomb interactions into ac- 
count within a tight-binding approach for BLG via a 
Hubbard model description of the carbon ir- electrons. In 
particular, since the neutrality point relates to half-filling 
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FIG. 1. (a) Bernal stacking of the honeycomb bilayer with 
intra- (inter)layer hopping t (t±) between the sublattices A, 
B and A', B' (A', B). Within the sublattices an equal number 
of sites have a coordination number z = 3 or 4. (b) Patching 
scheme of the Brillouin zone in the fRG. Dots denote the 
momenta at which the vertex function is evaluated. 



in the Hubbard model description, we take the opportu- 
nity to explore possible electronic instabilities using un- 
biased quantum Monte Carlo (QMC) methods. Our sim- 
ulations are furthermore augmented by functional renor- 
malization group (fRG) calculations 
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allows us to investigate the stability of the AF state ob- 
tained with QMC simulations over a broad range of the 
interaction strength. We find that within the AF ground 
state a local spin moment's participation in the AF order 
anticorrelates to its lattice coordination number z, with 
z = 3 or 4 for the Bernal stacking, an effect that we show 
to hold over the full parameter range from weak to strong 
electronic correlations. 

In the following we consider the Hamiltonian H = H$+ 
Hi n t, with the local interaction term H lnt — U Hj^rij^ 
and m, a = c\ a c i a the density operator at site i for 
spin a. Furthermore, Hq denotes the free tight-binding 
model Q containing both intralayer nearest neighbor 
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hopping t as well as inter layer hopping t±, as illustrated 
in Fig. 1(a). For the onsite interaction U and also a 
finite set of nonlocal density-density interaction param- 
eters, ah initio calculations list values for graphene and 
graphite (2(| , which have been used to explore the phase 
diagram of the honeycomb bilayer by means of the fRG 
approach that is also employed in the present work [22| . 
It was found that AF order is the dominant instability for 
interaction parameters with a shorter range than those 
for single layer graphene. As for a Bernal-stacked bilayer, 
where screening is expected to be effective, the antifcr- 
romagnet seems to be a viable candidate for the ground 
state of BLG. In this Letter, we employed an improved 
fRG patching of the Brillouin zone into 48 sectors, cf. 
Fig. 1(b), in order to obtain more accurate estimates for 
the critical energy scale A c where an electronic instabil- 
ity emerges during the fRG flow while successively inte- 
grating out the high-energy modes. Details on the fRG 
approach have been presented in Ref . I22I. From analyzing 
the structure of the resulting interaction vertex, we can 
identify the leading instability below A c for varying sets 
of initial coupling parameters pjj . For a broad range of 
pure Hubbard interactions U, we observe a flow to strong 
coupling with the signature of an AF instability and an 
exponential dependence of A c on U as discussed below. 
Up to an order of magnitude, A c can serve as estimate 
for the single particle gap A sp in the AF state. The AF 
instability is robust with respect to variations of the band 
structure, in particular the interlayer coupling, which we 
have explicitly checked for t± = O.li and t± = t (for 
BLG t± w 0.13i [III]). We take this as further motiva- 
tion for a systematic analysis of the local Hubbard model 
at half-filling. 

Our main findings result from analyzing this local Hub- 
bard limit, where we can efficiently employ a projec- 
tor QMC approach [2^], to perform a numerically ex- 
act evaluation of the ground state properties. We fur- 
thermore fix t± = t; while this takes us beyond the 
regime of realistic parameters for BLG, the choice for 
ij_ allows us to reliably study electronic instabilities, 
due to the well pronounced quadratic band touching 
at the Fermi level [2^|. We performed QMC simula- 
tions on finite systems of linear extent L (the number 
of sites being ./V = 4L 2 ) for L up to 12 with periodic 
boundary conditions. Our implementation also allows 
for the efficient measurement of unequal-time correlation 
functions (30| . From a fit of the imaginary-time dis- 
placed Green's function G(q, r) = (± J2 s ,a c \sa( T ) c ^a) 
to its long-time behavior, lim^oo G(q, r) cx e -TAsp(<5 > , 
the single-particle gap A sp = A sp (K) can be extracted 
without the need of an analytical continuation. Here, 
s labels the four orbitals per unit cell of the honey- 
comb bilayer. In order to gain information on the in- 
fluence of finite size (FS) effects we found it useful to 
compare also to fRG as well as to a AF mean-field the- 
ory (MFT) decoupling of the interaction term H lnt - We 
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FIG. 2. Single particle gap A sp from QMC simulations for 
different system sizes and the finite size extrapolation to the 
TDL using a polynomial fit function, along with the fRG crit- 
ical scale A c as a function of the local Coulomb repulsion U/t. 
The inset on the left shows the same data vs t/U in a semilog 
scale, exhibiting an exponential onset of the gap in the large 
t/U range. The inset on the right shows the fRG data along 
with MFT results for system sizes L = 129, 258, and 516. 

solved the resulting saddle-point equations self consis- 
tently on finite lattices and zero temperature to obtain 
the MFT order parameters for the sublattice magnetiza- 
tion and the associated Hartree-Fock single particle gap 
A™ FT = min {k} ^e(k) 2 + {Urn) 2 = Urn, where e(k) the 
single particle dispersion of the free Hamiltonian Hq. 

In Fig. [2j we present our results for the single parti- 
cle gap A sp obtained from QMC simulations and MFT 
alongside the fRG critical scale A c as a function of U/t. 
The QMC values exhibit a continuously increasing A sp 
for all system sizes. The FS extrapolation to the ther- 
modynamic limit (TDL) using a second order polynomial 
yields a continuous onset in the thermodynamic limit. Fi- 
nite size extrapolation of the available data points decep- 
tively suggest a finite critical value of U for the transition 
from the semimetal to the Mott insulator. We will show 
in the following that this can be attributed to pronounced 
FS effects at low energies. The fRG critical scale A c (open 
circles) for the same parameters indeed reproduces such 
a continuously increasing associated single particle gap, 
for all finite values of U. 

To identify whether the observed gap indeed is the con- 
sequence of a U = + instability, we plot A sp as a func- 
tion oit/U in the insets of Fig. [5] The fRG data show an 
exponential opening of A sp in the large t/U regime (left 
inset). Accordingly, the FS extrapolated QMC data fol- 
low the same behavior. One can readily see that larger 
lattices are needed in order to clearly identify this ex- 
ponential onset at smaller values of U/t. The same ef- 
fect is observed already within the MF approach (right 
inset): increasingly larger system sizes allow us to iden- 
tify the exponential opening of the single particle gap 
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FIG. 3. The sublattice magnetization vs. U/t for sites with 
coordination number z = 3 (open symbols) and z = 4 (filled 
symbols) from (a) MFT in the TDL, (b) QMC simulations for 
different system sizes, and (c) in the Heisenberg limit vs. the 
interlayer AF exchange coupling J± / J. The data in (c) result 
from SSE QMC simulations, after extrapolation to the TDL; 
also shown in (c) is the overall staggered magnetization m. 
Part (d) shows the local dynamical spin susceptibility Xo-(i^) 
for sites with z = 3 and 4 for L = 12, U/t = 4. 



A sp oc exp(— at/U). From our fRG data we can extract 
the exponent a « 16, which is close to the Hartree-Fock 
value of 9nt 2 /2tj_ [31|. Deviations from this exponen- 
tial behavior emerge beyond an intermediate coupling 
strength of U/t « 2 and relate to the onset of the strongly 
correlated regime. Note that within this model, the en- 
ergy gap for realistic values of the Hubbard U ~ 3t are 
rather large, A sp « 0.07t « 200 mcV 

We identify the order parameter associated with the 
single particle gap from the low-energy vertex in the fRG 
and from corresponding correlation functions in QMC 
simulations, consistently, to be long-range AF order, cor- 
related between both layers. To quantify this order 
within QMC simulations, we measure the overall stag- 
gered structure factor Saf = J2i j e i e j(^i ' Sj) from 
which we obtain the mean staggered magnetization per 
lattice site as m = y 7 Saf/ (4£ 2 )- Here, = ±1 if site 
i belongs to the magnetic sublattice A, A' (B, B'), as 
indicated by the white (black) spheres in Fig. 1. In or- 
der to probe in more detail the magnetic correlations, 
we also consider the following restricted structure factors 
for sites with coordination numbers z = 3 and 4 (for a 
system of linear size L, there are 2L 2 such sites each): 



tain local order parameters m z = yl Saf,z/ (2£ 2 ) for lat- 
tices sites with z = 3,4. While the overall staggered 
magnetization m steadily increases with U > like the 
single particle gap, we observe pronounced differences in 
the two sublattice magnetizations, which arise due to the 
presence of inequivalcnt sites in the lattice structure [cf. 
Fig. EJa)-(b)]. In particular, we find that sites with the 
higher coordination z = 4 (filled symbols) exhibit a lower 
ordered moment, at odds with the usual intuition that 
high coordination favors more robust Neel order. An in- 
crease of z on one sublattice will generically make the or- 
dering more robust everywhere, although not uniformly 
so for all sublattices. The same hierarchy of magnetic 
moments can also be inferred from fRG calculations for 
a wider range of nonlocal interactions by comparing the 
relative strengths of the effective interactions on the dif- 
ferent sites of a unit cell. 

Similar effects of an enhanced magnetic order near low- 
coordinated sites have previously been observed in local- 
ized quantum s pin models on other inhomogeneous lat- 
tice structures [32J. Here, the joint bonds on sites with 
coordinated number z = 4 interconnect the two layers [cf. 
Fig. [1] (a)]. With increasing interaction U and hopping 
t±, the moments along these bonds hybridize between 
the layers and tend to form spin singlets, suppressing the 
participation in the long range AF order on these sites. 
In Fig. 0(c) , wc consider the staggered magnetizations 
in the large-C/ limit, wherein the model becomes a spin- 
only Heisenberg model, with an intralaycr exchange cou- 
pling J = At 2 /U and an interlayer coupling J± = At\/U . 
To study this Heisenberg limit of the Hubbard model, 
we employed the stochastic series expansion (SSE) QMC 
approach [3^, |34[ i and present our results after an extrap- 
olation to the TDL. We find that all three order param- 
eters exhibit an initial increase upon increasing Ji/J. 
Furthermore, while 7713 saturates for large J±/J, m and 
7TI4 scale to zero in the large J± limit. These results show 
that for all finite values of J/J±, the system remains an- 
tifcrromagnctically ordered, but with a suppressed stag- 
gered moment on the z — 4 sites for large Jj_/J, due to 
the strong tendency towards forming J±- singlets along 
these interlayer bonds. If the two honeycomb lattices 
were stacked such that each lattice site would be coupled 
via J± to the other layer, a complete decoupling into lo- 
cal singlets would destroy the AF order beyond a finite 
critical value of J±/J. 

Returning to the Hubbard model, the local dynamical 
spin susceptibility x<r,i(u) = J2n IH 5 ? |0}| 2 S(E n - E - 
ui) also exhibits defined differences depending on the lo- 
cal coordination, cf. the QMC data in Fig. [3jd). While 
the peak near w = 0, related to the Goldstone mode in 
the TDL, is shared by both types of sites, the z = 4 sites 
exhibit a considerable shift of the residual spectral weight 
to larger energies as compared to the z = 3 sites, in ac- 
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FIG. 4. Finite size behavior of (a) the spin gap A CT and (b) the 
AF structure factor Saf,3 on lattice sites with coordination 
number z = 3. The inset in (b) illustrates the discretization 
effects on the free dispersion relation, a dominant source of 
finite size effects. 



cordance with the associated reduced magnetic moment 
for z = 4. 

To explore the global spin dynamical properties, 
Fig. Eta) shows the spin gap A CT from QMC simulations 
as a function of U/t for different system sizes (we obtain 
A CT from the time-displaced spin-spin correlation func- 
tion in the AF sector, S A f(t) = jj J^i.j e i e j(Si( T ) ' s i))- 
The data for finite lattices exhibits a pronounced peak at 
intermediate values of U/t; for increasingly larger lattice 
sizes, the position of this spin-gap dome shifts towards 
lower values of U/t, and its magnitude decreases, sug- 
gesting that in the TDL the spin gap vanishes for all val- 
ues of U/t. While this is consistent with the emergence 
of Goldstonc modes which originate in the spontaneous 
breaking of the SU(2) spin symmetry in the AF phase, it 
furthermore illustrates the pronounced FS effects on the 
accessible range of system sizes. Similarly distinct FS 
corrections are also evident in the AF structure factor 
5af,3, shown vs. 1/L in Fig. EJb). Consider, for exam- 
ple, the data for U/t = 2.4: the initial downscaling of 
Saf,3 on small lattice sizes would suggest a magnetically 
disordered state. However, at larger lattice sizes the scal- 
ing behavior changes, and eventually a strong increase of 
Saf,3 with system size accounts for the formation of long- 
range AF order in the TDL. This peculiar FS scaling in 
fact arises both in the QMC simulations and the MFT 
order parameters (not shown). The FS effects become 
more pronounced for smaller t±. The data in Fig. E£b) 
also show that the corresponding crossover length scale 
beyond which the eventual increase of Sap, 3 sets in, in- 
creases with decreasing values of U ft. Finite lattices, i.e., 
momentum space discretization, introduce a correspond- 
ing artificial gap which acts as a cutoff for correlations — 
an issue which afflicts all finite lattice simulations. When 



dealing with a Fermi surface instability as in the present 
case, this generic FS effect becomes predominant on the 
accessible system sizes. The inset of Fig. |4|b) illustrates 
this discretization of the free dispersion around the Fermi 
level. Only for sufficiently large lattices will the parabolic 
band be approximated to such an extent as to probe the 
TDL nonlinear low energy dispersion. 

In conclusion, we found from quantum Monte Carlo 
simulations, that the Hubbard model on the Bernal- 
stacked honeycomb bilaycr as a basic model for BLG is 
prone to a Fermi point instability, which triggers layered 
AF order. Characterized by their different coordination 
numbers, sublattice sites sustain different magnetic mo- 
ments. This peculiar local structure of the AF state re- 
lates to its stability in the strong interlayer tunneling 
region. A full quenching of the magnetic moments on 
the z = 4 sites emerges only in the (unrealistic) strong 
interlayer coupling limit and would eventually realize the 
layered AF state observed within chiral two-band models 
for bilayer graphenc. Functional renormalization group 
calculations support the stability of the AF state and its 
moment distribution over a wide range of coupling pa- 
rameters. In case the experimental support for a gapped 
state in BLG will be substantiated, it might be inter- 
esting to search for such an inhomogeneous AF state by 
local moment-sensitive probes such as magnetic scanning 
tunneling microscopy. 
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